/*********************************************************
 Main Window! 
 
 Click the "Run" button to Run the simulation.
 
 Change the geometry, flow conditions, numercial parameters
 visualizations and measurments from this window. 
 
 *********************************************************/

BDIM flow;
//BodyUnion3 body;
CircleBody body;
FloodPlot flood;
MovieMaker mm;
SaveData stream;
SaveData transverse;
SaveData vorticity;
SaveData vortGrad;
SaveData drag;
boolean recording = true;
int n=(int)pow(2, 10)+2; // number of grid points along a side
float t=0, D = n/10, dtm=0;
float dR = .125;

void setup() {
  float xc = n/2.5, yc=n/2;

  float theta = 60*3.1415926/180;

  float gR = .05;
  float r=(D/2+gR*D+dR*D/2);

  // set window view area and zoom
  int area = 300000;
  float s = sqrt(area);
  size(int(s), int(s));
  float zoom = 3;
  Window view = new Window(int((n-n/zoom)/2), int((n-n/zoom)/2), int(n/zoom), int(n/zoom));

  // create circle with rotating cylinders
//  body = new CircleBody(xc, yc, D, view);
//  body = new BodyUnion3( new CircleBody(xc, yc, D, view), 
  body = new CircleBody(xc+r*cos(theta), yc+r*sin(theta), dR*D, view);
//  new CircleBody(xc+r*cos(theta), yc-r*sin(theta), dR*D, view));

  float Re = 500;
  float nu = D/Re;
  flow = new BDIM(n, n, 0, body, nu, true);   // QUICK with adaptive dt

  // create vorticity heat map
  int nlines = int(30*zoom);        // number of streamlines
  flood = new FloodPlot(view);
  flood.range = new Scale(-1, 1);
  flood.setLegend("vorticity");

// to resume from a save file, put the file in the "data" folder
//  resumeBDIM("afterTransients_xi0_res10.dat");

// create data files
//  stream = new SaveData("stream.txt", body.coords, 1, 1, 1, 1);
//  transverse = new SaveData("transverse.txt", body.coords, 1, 1, 1, 1);
//  vorticity = new SaveData("vorticity.txt", body.coords, 1, 1, 1, 1);
//  vortGrad = new SaveData("vortGrad.txt", body.coords, 1, 1, 1, 1);
  drag = new SaveData("drag_control_Re500.txt", body.coords, int(D/nu), 1, 1, 1);
}


void draw() {
  t += flow.dt/D;

  //  if (t>125) {
  //    body.c.rotate(2*flow.dt*flow.dt*2/(dR*D));
  //    body.b.rotate(-2*flow.dt*flow.dt*2/(dR*D));
  //  }

  body.update();
  flow.update(body);
  flow.update2();
  flood.display(flow.u.vorticity());
  body.display();
  flood.displayTime(t);

  drag.addSimple(t, body.pressForce(flow.p).x*2/D);

//  if (recording) {
//    if (t>60) {
//      if (t>dtm) {
//        stream.addField(t, n, flow.u.x);
//        transverse.addField(t, n, flow.u.y);
//        dtm = t+.1;
//      }
//    }
//  }

  if (t>40) {
//    stream.finish();
//    transverse.finish();
//    vorticity.addField(t, n, flow.u.vorticity());
//    vortGrad.addField(t, n, flow.u.vorticity().gradient().y);
//    vortGrad.finish();
//    vorticity.finish();
    drag.finish();

    // save flow state
    writeBDIM("afterTransients_xi0_res10.dat");

    exit();
  }
}


void writeBDIM( String name ) {
  PrintWriter output;
  output = createWriter(name);
  output.println(t);
  output.println(flow.dt);
  for ( int i=0; i<n; i++ ) {
    for ( int j=0; j<n; j++ ) {
      output.println(""+flow.u.x.a[i][j]+", "+flow.u.y.a[i][j]+", "+flow.p.a[i][j]);
    }
  }
  output.flush(); 
  output.close();
}

void resumeBDIM( String name ) {
  float[] data;
  String[] stuff = loadStrings(name);
  t = float(stuff[0]);
  flow.dt = float(stuff[1]);
  for ( int i=0; i<n; i++ ) {
    for ( int j=0; j<n; j++ ) {
      data = float(split(stuff[2+i*n+j], ','));
      flow.u.x.a[i][j] = data[0];
      flow.u.y.a[i][j] = data[1];
      flow.p.a[i][j] = data[2];
    }
  }
}


void writeBDIM( String name ) {
  PrintWriter output;
  output = createWriter(name);
  output.println(t);
  output.println(flow.dt);
  for ( int i=0; i<n; i++ ) {
    for ( int j=0; j<n; j++ ) {
      output.println(""+flow.u.x.a[i][j]+", "+flow.u.y.a[i][j]+", "+flow.p.a[i][j]);
    }
  }
  output.flush(); 
  output.close();
}

void resumeBDIM( String name ) {
  float[] data;
  String[] stuff = loadStrings(name);
  t = float(stuff[0]);
  flow.dt = float(stuff[1]);
  for ( int i=0; i<n; i++ ) {
    for ( int j=0; j<n; j++ ) {
      data = float(split(stuff[2+i*n+j], ','));
      flow.u.x.a[i][j] = data[0];
      flow.u.y.a[i][j] = data[1];
      flow.p.a[i][j] = data[2];
    }
  }
}

